data <- clean_LISS_data

# Adjust Variable Types & Add Unit Weights
data <- data %>%  mutate(high_flex = as.numeric(high_flex), high_telecommute = as.numeric(high_telecommute), high_meaning = as.numeric(high_meaning))
data <- data %>% mutate(weight = 1)

# Select Fulltime
data <- subset(data, whcat == 2 | whcat == 1)

# Create Subdatasets - Men & Women
data_mothers <- subset(data, woman == 1)
data_fathers <- subset(data, woman == 0)

# Function for weighted mean with CI
source(file.path(dir, "DataCode/Code/Functions/weighted_mean_ci.R"))

# Amenities by Gender
mean_flex <- weighted_mean_ci(data, high_flex, weight, woman)
mean_tel <- weighted_mean_ci(data, high_telecommute, weight, woman)
mean_meaning <- weighted_mean_ci(data, high_meaning, weight, woman)

mean_summary <- mean_flex %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_tel %>% rename(tel_mean = mean, tel_lb = lb, tel_ub = ub), by = "woman") %>%
  inner_join(mean_meaning %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "woman")

# Amenities by children: Mothers
mean_flex_m <- weighted_mean_ci(data_mothers, high_flex, weight, children)
mean_tel_m <- weighted_mean_ci(data_mothers, high_telecommute, weight, children)
mean_meaning_m <- weighted_mean_ci(data_mothers, high_meaning, weight, children)

mean_summary_m <- mean_flex_m %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_tel_m %>% rename(tel_mean = mean, tel_lb = lb, tel_ub = ub), by = "children") %>%
  inner_join(mean_meaning_m %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "children")

# Amenities by children: Fathers
mean_flex_f <- weighted_mean_ci(data_fathers, high_flex, weight, children)
mean_tel_f <- weighted_mean_ci(data_fathers, high_telecommute, weight, children)
mean_meaning_f <- weighted_mean_ci(data_fathers, high_meaning, weight, children)

mean_summary_f <- mean_flex_f %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_tel_f %>% rename(tel_mean = mean, tel_lb = lb, tel_ub = ub), by = "children") %>%
  inner_join(mean_meaning_f %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "children")

# Prepare for Plots
plot_data <- mean_summary %>%
  pivot_longer(cols = c(flex_mean, tel_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(woman == 1, 0.10, 0),
      variable == "flex_mean" ~ if_else(woman == 1, 1, 0.90),
      variable == "tel_mean" ~ if_else(woman == 1, 1.90, 1.80)
    ),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb,
      variable == "tel_mean" ~ tel_lb),
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub,
      variable == "tel_mean" ~ tel_ub)
  )

plot_data_m <- mean_summary_m %>%
  pivot_longer(cols = c(flex_mean, tel_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(children == 1, 0.35, 0.25),
      variable == "flex_mean" ~ if_else(children == 1, 1.25, 1.15),
      variable == "tel_mean" ~ if_else(children == 1, 2.15, 2.05)
    ),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb,
      variable == "tel_mean" ~ tel_lb,
    ),
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub,
      variable == "tel_mean" ~ tel_ub)
  )

plot_data_f <- mean_summary_f %>%
  pivot_longer(cols = c(flex_mean, tel_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(children == 1, 0.60, 0.50),
      variable == "flex_mean" ~ if_else(children == 1, 1.50, 1.40),
      variable == "tel_mean" ~ if_else(children == 1, 2.40, 2.30)
    ),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb,
      variable == "tel_mean" ~ tel_lb),
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub,
      variable == "tel_mean" ~ tel_ub
    )
  )

# Color Scheme
bw_palette <- c("grey10", "grey30", "grey50", "grey70", "grey85", "grey95")

# Sort for plot
plot_data = arrange(plot_data, by_group = amenity)
plot_data_f = arrange(plot_data_f, by_group = amenity)
plot_data_m = arrange(plot_data_m, by_group = amenity)


plot_amenities_LISS <- ggplot() + theme_bw() +
  geom_col(data = plot_data, aes(x = amenity, y = mean, fill = rep(c("Men", "Women"),3)), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  geom_col(data = plot_data_m, aes(x = amenity, y = mean, fill = rep(c("No Mother", "Mother"),3)), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data_m, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  geom_col(data = plot_data_f, aes(x = amenity, y = mean, fill = rep(c("No Father", "Father"),3)), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data_f, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  xlab("") + ylab("Percentage with Amenity") +
  scale_x_continuous(breaks = c(0.30, 1.25, 2.20), labels = c("Work Meaning", "Schedule Adaptability", "Telecommuting")) +
  theme(
    legend.key = element_rect(fill = "white", colour = "white"), 
    legend.position =  "bottom", 
    legend.box = "horizontal", 
    legend.text = element_text(size = 12),
    axis.text.x = element_text(size = 14, color = "black"),  # Increase x-axis label size
    axis.title.y = element_text(size = 14),  # Increase x-axis label size
    axis.text.y = element_text(size = 14, color = "black"),  # Increase y-axis label size
    axis.ticks.x = element_blank()  # Remove x-axis ticks
  ) +
  labs(color = NULL, fill = NULL) +  # Set the fill legend title to NULL
  scale_fill_manual(values = bw_palette, breaks = c("Men", "Women", "No Mother", "Mother", "No Father", "Father"))  + 
  guides(color = guide_legend(override.aes = list(size = 7), ncol = 6, name = NULL))

save_PlotAmenities_LISS <- file.path(graph_dir, "AmenityHeterogeneityGenderChild_LISS_FT.pdf")

# Save as PNG
pdf(save_PlotAmenities_LISS) #, width = 800, height = 600)
print(plot_amenities_LISS)
dev.off()

# Remove newly created dataframes
rm("data", "data_mothers", "data_fathers", "mean_flex", "mean_tel", "mean_meaning", "mean_summary", "mean_flex_m", "mean_tel_m", "mean_meaning_m",  "mean_summary_m", "mean_flex_f", "mean_tel_f", "mean_meaning_f",  "mean_summary_f")
